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Abstract 

Lattice models describing the spatial spread of rabies among foxes are studied. 
In these models, the fox population is divided into three-species: susceptible (S), 
infected or incubating (I), and infectious or rabid (R). They are based on the fact 
that susceptible and incubating foxes are territorial while rabid foxes have lost their 
sense of direction and move erratically. Two different models are investigated: a one- 
dimensional coupled-map lattice model, and a two-dimensional automata network 
model. Both models take into account the short-range character of the infection 
process and the diffusive motion of rabid foxes. Numerical simulations show how 



the spatial distribution of rabies, and the speed of propagation of the epizootic front 
depend upon the carrying capacity of the environment and diffusion of rabid foxes 
out of their territory 



1 Introduction 



Rabies is one of the oldest recorded infectious diseases. Records from the Middle Ages 
show that rabies was then widespread in western Europe. Throughout the world, dogs are 
important vectors, however, there exist many early references to several species of wildlife 
implicated in the spread of the disease U, @, ||. Nowadays, most cases of rabies are 
found among animals (dogs, cats, livestock, wildlife). According to Macdonald ||, the 
present European epizootic (epidemic among animals) began south of Gdansk (Poland) 
in 1939. Moving westwards, it was first recorded in France in 1968, the Netherlands in 
1974, Spain in 1975 and Italy in 1977. Its main vector and victim is the red fox. 

Rabies is an acute infectious disease of the central nervous system caused by a virus. 
The disease is transmitted from rabid to susceptible foxes, usually by biting. About half 
of the infectious foxes become agressive, and lose their sense of direction and territorial 
behavior, wandering randomly. 

Murray and co-workers 0, f|, f| studied different models of the spread of rabies 
among foxes formulated in terms of partial differential equations. In these models, the 
fox population is divided into three species: susceptible, infected but non-infectious, and 
rabid foxes which are infectious. The principal assumptions are: 

(i) Susceptibles evolve in time according to a logistic equation. 

(i) Susceptibles become infected at an average rate per capita proportional to the num- 
ber of rabid foxes present. 

(iii) Infected foxes become rabid after an average incubation period. 

(iv) Rabid foxes die after an average duration of the disease. 

(v) Rabid and infected foxes also die of causes other than rabies. 

(vi) The spatial spread of the disease is essentially due to the random motion of the 
rabid foxes 

Such models have unquestionably contributed to our understanding of the spread of an 
infectious disease, but they do not take correctly into account the short-range character of 



the infection process, and neglect spatial correlations. In order to include these features, 
we have built up and studied two lattice models in which the spread of the epidemic is 
viewed as the growth of a random cluster on a lattice. We will first describe a coupled-map 
lattice model, and then a two-dimensional automata network model. 

2 Lattice models 

The spread of an epidemic in a population is a complex process. When modeling such 
a process, among the many features which are likely to be important, we should, how- 
ever, only retain the few relevant ones which are thought to play an essential role in the 
interpretation of the observed phenomena. In both models, at each time step, the fox 
population evolves according to the following rules: 

(1) A susceptible has a probability b to give birth to a susceptible at a nearby empty 
location. Infected and rabid foxes do not give birth. 

(2) Susceptibles and infected have a probability d to die due to natural causes, and dg 
multiplied by the total local population density of neighboring foxes {%. e. susceptibles 
+ infected + infectious) due to lack of food. 

(3) Susceptibles become infected at a rate proportional to the local density of neighbor- 
ing rabid foxes, the proportionality factor is the probability to be infected. 

(4) Infected, i.e. incubating, foxes become rabid with a probability p r . 

(5) Non-rabid foxes being territorial, susceptibles and infected foxes evolve without 
moving to a neighboring territory. 

(6) Rabid foxes have a probability d r to die due to the disease, and dg multiplied by the 
total local density of neighboring foxes due to lack of food. 

(7) Rabid foxes move at random to a neighboring location. 

In order to exhibit some realistic features of the spread of a rabies epidemic, follow- 
ing the analysis of Murray et al B, our models have to contain a minimum number of 



parameters. There is one source term coming from the birth of susceptibles (parameter 
b). To account for the three death processes, natural causes, lack of food, and rabies, we 
need 3 parameters (d, dg, and d r ). We need to introduce a parameter which measures 
the probability to be infected by contact (parameter p^. The existence of an incubation 
period is an essential feature to exhibit decreasing periodic fluctuations following the main 
wave front of the susceptible population (parameter p r ). 

2.1 Coup led- map lattice model 

A coupled-map lattice is a dynamical system in which space and time are discrete variables 
while states are continuous J7|, ||. Here the state of the system is represented by the 
function P : (i,t) i— > P(i,t), where (i,t) 6ZxN, and P G [0, l] 3 is a three-dimensional 
vector whose components are S(i,t), I(i,t), and R(i,t) which denote, respectively, the 
densities at site i and time t of susceptible, infected and rabid foxes. Z is the set of all 
integers and N is the set of nonnegative integers. In this model, each site corresponds to a 
specific territory. Susceptibles and infected foxes evolve without moving to a neighboring 
site. On the contrary, rabid foxes, which have lost their sense of direction, move to one of 
their two neighboring sites with equal probabilities . According to our assumptions (Rules 
1 to 7), the dynamics of the system is governed by the following recurrence relations: 

S(i, t + l) = S(i, t) - dS(i, t) - d e N(i, t)S(i, t) 

+ 6(1 - N(i, t))S(i, t) - Pt R(i, t)S(i, t) 
t + 1) = t) - dl(i, t) - d e N(i, t)I(i, t) 
+ PiR(i,t)S(i,t) -p r I(i,t) 

R(i, t + l) = R(i, t) - d r R(i, t) - dtN(i, t)R(i, t) + p r I(i, t) 
+ D(R(i -l,t)) + R{i + 1, t) - 2R{i, t)), 

where N(i, t) = S(i, t) + t) + R(i, t) is the total fox density at site i and time t, and 
D is the diffusion coefficient characterizing the random motion of the rabid foxes. 



2.2 Automata network model 

An automata network is a fully discrete dynamical system. That is, space, time and states 
are discrete variables. In our model, the space consists of a square L x L. Since we are 
interested in the spread of an epidemic starting a the center of the lattice, we did not 
choose cyclic boundary conditions. As for a coupled-map lattice, the time variable is a 
nonnegative integer. The state of the system at time t is described by the function s t : 
i— > s(i,j; t), where E 1? L , t EN, and s(i, j; t) can take four values corresponding 
to the four possible states of the site since this site can be either empty or occupied 
by one of the three fox species (susceptible, infected, rabid). At each time step, the state 
of the system evolves according to the successive application of the two following subrules: 

• first a local rule describing birth, death and infection processes, which is applied to 
all sites synchronously; 

• then, a motion rule mimicking rabid fox erratic motion. 

Sub rule (1) is a probabilistic, two-dimensional, four-state cellular automaton rule, 
which is defined as follows. At each time step, 

(i) each susceptible has a probability b to give birth to a susceptible at an empty first- 
neighbor site, therefore, the probability that an empty site, having z s susceptibles 
among its four first-neighbors, becomes occupied by a susceptible is 1 — (1 — b) Zs ; 

(ii) each susceptible has a probability d to die of natural causes; 

(iii) each susceptible having Z; L rabid foxes among its four first-neighbors becomes infected 
with a probability 1 — (1 — Pi) Zl ; 

(iv) each infected has a probability d of die of natural causes; 

(v) each infected has a probability p r to become rabid; 

(vi) a rabid fox has a probability d r to die of the disease; 

(vii) each fox (either susceptible, or infected, or rabid) having Zf foxes among its first- 
neighbors has a probability dgZf to die due to lack of food. 



Subrule (2) can be described as follows: At time t, a rabid fox is selected at random 
to perform a move to a first-neighbor site also selected at random. If the site is empty 
the fox will effectively move otherwise it will not. This process is repeated miV^(£)L 2 , 
where Nn(t) denotes the density of rabid foxes at time t. The parameter m represents 
the average number of tentative move per rabid fox at time t. Note that this sequential 
diffusive process allows some rabid foxes to move more than others. 

Automata networks of this type are called diffusive cellular automata. Their general 
properties have been studied by Boccara et al |§. They have been used to build up 
various epidemic models jHj in which the motion of the individuals play an important 
role in the spread of the epidemic. 



3 Numerical simulations: results and analysis 
3.1 Coup led- map lattice model 

In our simulations, the lattice size is 441. The sites are labeled from —220 to 220. The 
initial susceptible fox densities have the same value at all sites i, while the initial infected 
and rabid fox densities are nonzero only at the central site. These values are: S(i, 0) = 0.6. 
1(0, 0) = 0.005, and R(0, 0) = 0.005. For all t <E N, at the boundaries of the chain, the 
densities of the various fox species satisfy the condition 

S(L + 1, t)+I(L + l,t) + R(L + l,t) = 

S(-L - 1, t) + I(-L - 1, t) + R(-L - 1, t) = 

As reported by Macdonald [|3| and Murray ||, the epidemic spread of rabies among 
foxes is characterized by a traveling epizootic wave front followed by periodic decreas- 
ing fluctuations of susceptibles density which tends to its steady state. Our numerical 
simulations show that our coupled-map lattice model clearly exhibits these features, as 
illustrated in Figure 1 which represents, at a given time, the variations of susceptible, and 
rabid fox densities, as a fonction of site location for two different values of dg. Increas- 
ing di decreases the height of the peak of the first and subsequent outbreaks of infected 



and rabid foxes. A smaller value of dg is equivalent to a larger carrying capacity of the 
environment, therefore, a larger carrying capacity implies a more severe epidemic. 

Figure 2 illustrates the influence of the diffusion coefficient D. As expected, increasing 
D increases the speed at which the epidemic spreads (see also Figure 3). While the values 
of the amplitudes of the various fox densities do not change, the distance between two 
successive outbreaks increases with D. 

According to epidemiological evidence (cf. references in Murray ||) there exists a 
threshold value for the carrying capacity below which rabies die out. In the case of our 
model we found that for dg > 0.0456 the epizootic does not spread (see Figure 4). This 
threshold depends, of course, upon the values of the other parameters. 

We have determined numerically the speed of the epizootic wavefront as a function of 
the diffusion coefficient D for two different values of the parameter dg. Our results are 
represented in Figure 3. As expected from a dimensional argument, this speed varies as 
\f~D. Increasing dg, we verify again that the speed of the epizootic wavefront decreases. 

We have also determined numerically the speed of the epizootic wavefront as a function 
of the parameter dg. Figure 4 shows that, as we already mentioned, this speed is a 
decreasing function of dg, which goes to zero at a threshold value. 

3.2 Automata network model 

In our simulations, the lattice size is 201 x 201. Simulations start from a random initial 
configuration. The initial densities of susceptible, infected and rabid foxes are, respec- 
tively, equal to 0.6, 0.005, and 0.005. Rabid and infected foxes exist only inside a disk 
of radius 10 in the initial configuration. Our results are averages over 50 to 100 different 
initial configurations. 

As for the coupled-map lattice model, we have studied the influence of the carrying 
capacity and the diffusion on the various fox species. We found similar results as illustrated 
in Figures 5a-5c. Increasing dg decreases the height of the peak of the first and subsequent 
outbreaks of rabid foxes (Figures 5a and 5b), while increasing m increases the speed at 
which the epidemic spreads (Figures 5b and 5c). Above a threshold value of dg, which 
is equal to 0.053 for m = 0.5 and 0.0245 for m = 0.1, the epizootic does not spread (see 



Figure 9). 

As time increases, rabid foxes can be found at larger distances from the center of the 
lattice. For different values of the parameter m, we have determined the fractal dimension 
of the rabid cluster — that is the cluster containing all rabid foxes — as a function of time. 
If W is the radius of this cluster, its fractal dimension 5 is defined by N R (t) = vrr^ ax , 
where N R (t) is the number of rabid foxes in the cluster. Figure 6 shows that 5 tends to 
a constant value as time increases. This limit value increases when m decreases since, as 
expected, the density of rabid foxes if higher for lower values of m. For de respectively 
equal to 0.005, 0.01, and 0.015, the corresponding limit values are 1.27, 1.25, and 1.22. 

Figures 7a- 7f show how the cluster containing the infected and infectious foxes grows as 
a function of time. Infected and infectious foxes are essentially localized at the boundary 
of the cluster. This feature corresponds to the small subsequent outbreaks following the 
first one. This is also clearly illustrated in Figures 5. 

As for the coupled-map lattice model, we have determined numerically how the speed 
of the epizootic wavefront varies with the carrying capacity and the diffusion. Our results 
are represented in Figures 8 and 9. The curves representing the variation of the speed of 
the epizootic wavefront as a function of de show the existence of a threshold value, which, 
in particular, depends upon the parameter m. 

4 Conclusion 

We have investigated two different models of the spatial spread of rabies among foxes: 
a one-dimensional coupled-map lattice model, and a two-dimensional automata network 
model. In both models, the fox population is divided into three-species: susceptible, in- 
fected or incubating, and infectious or rabid. They are based on the fact that susceptible 
and incubating foxes are territorial while rabid foxes have lost their sense of direction 
and move erratically out of their territory propagating the disease, coupled-map lattice 
models and automata networks models have the advantage, compared to models formu- 
lated in terms of differential equations, to take into account the short-range character of 
the infection process. We have essentially studied how the spatial distribution of rabies, 



and the speed of propagation of the epizootic front depend upon the carrying capacity 
or food availability of the environment and parameters characterizing the erratic motion 
of rabid foxes out of their territory. In agreement with ecological studies, our numerical 
simulations show that, decreasing food availability slows down the spread of the disease, 
and, that below a certain threshold, rabies eventually dies out. On the other hand in- 
creasing the parameters measuring the diffusive motion of rabid foxes favors the spread 
of the disease. 



References 



[1] G. M. Baer, The natural history of rabies, (New York: Academic Press, 1975). 

[2] R. K. Sikes Sr, Rabies in infectious diseases of wild animals, eds. J. W. Davis, L. 
H. Karstad and D. O. Trainer, (Ames: The Iowa State University Press, 1981). 

[3] D. W. Macdonald, Rabies and wildlife: a biologist's perspective, (Oxford: Oxford 
University Press, 1980). 

[4] A. Kallen, P. Arcuri and J. D. Murray, J. Theor. Biol. 116, 377 (1985). 

[5] J. D. Murray, E. A. Stanley and D. L. Brown, Proc. Roy. Soc. B. 229, 111 (1986). 

[6] J. D. Murray, Mathematical biology, (Heidelberg: Springer- Verlag, 1989). 

[7] K. Kaneko, Prog. Theo. Phys. 72, 480 (1984); K. Kaneko, Prog. Theo. Phys. 74, 
1033 (1985), K. Kaneko, Chaos, 2, 279 (1993). 

[8] I. Waller and R. Kapral, Phys. Rev. A 30, 2047 (1984), R. Kapral, Chemical waves 
and coupled-map lattices, in Theory and applications of coupled-map lattices, ed. K. 
Kaneko, (New York: Wiley, 1993) 

[9] N. Boccara and M. Roger, Int. J. Mod. Phys. C 5, (3) 581 (1994); N. Boccara, J. 
Nasser and M.Roger, Int. J. Mod. Phys. C 5, (3) 537 (1994). 

[10] N. Boccara and K. Cheong, J. Phys. A 25, 2447 (1992); N. Boccara and K. Cheong, 
J. Phys. A 26, 3707 (1993); N. Boccara, K. Cheong and M. Oram, J. Phys. A 
27, 1585 (1993); N. Boccara, Automata networks models of interacting populations, 
in Cellular automata, dynamical systems, neural networks, eds E. Goles and S. 
Martinez, (Dordrecht: Kluwer, 1994). 



Figure captions 

Figure 1- Susceptible (S) and rabid (R) densities as functions of site location at time 
t = 1800 for b = 0.01, d = 0.001, d r = 0.1, pi = 0.8, p r = 0.05, D = 0.1, d e = 0.01 
(continuous line), and de = 0.015 (dashed line). 

Figure 2- Fox populations densities as a function of site location at time t = 1800 for 
b = 0.01, d = 0.001, d e = 0.01, d r = 0.1, p x = 0.8, p r = 0.05, D = 0.1 (dashed line), and 
D = 0.2 (continuous line). 

Figure 3- Speed of the epizootic wavefront as a function of D for b = 0.01, d = 0.001, 
d r = 0.1, pi = 0.8, p r = 0.05, and two different values of de indicated in the figure. These 
two set of data can be fitted with functions of the form A\f~D (continuous lines) 
Figure 4- Speed of the epizootic wavefront as a function of de for b = 0.01, d = 0.001, 
d r = 0.1, p x = 0.8, p r = 0.05, and D = 0.1. 

Figure 5- Fox populations densities as a function of site location at time t = 1600 for 
b = 0.01, d = 0.001, d r = 0.1, pi = 0.8, p r = 0.05, and (a) d e = 0.01, m = 0.1, (b) 
d e = 0.015, m = 0.1, (c) de = 0.015, m = 0.2. 

Figure 6- Fractal dimension of the rabid cluster as a function of time for b = 0.01, 
d = 0.001, d T = 0.1, pi = 0.8, p r = 0.05, m = 0.2 and different values of de indicated in 
the figure. 

Figures 7- Spatial patterns obtained after different numbers of time steps. The symbols: 
' , A, and * denote, respectively, susceptible, infected and rabid foxes. Parameters values: 
b = 0.01, d = 0.001, d t = 0.001, d r = 0.1, p t = 0.8, p r = 0.05, m = 0.2. (a) t = 0, (b) 
t = 50 (c)t= 100, (d) t = 200, (e) t = 500, (f) t = 800. 

Figure 8- Speed of the epizootic wavefront as a function of m for b = 0.01, d = 0.001, 
d r = 0.1, pi = 0.8, p r = 0.05, and two different values of de indicated in the figure. 
Figure 9- Speed of the epizootic wavefront as a function of de for b = 0.01, d = 0.001, 
d r = 0.1, pi = 0.8, p r = 0.05, and two different values of m indicated in the figure. 
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